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Abstract 

The quantum many-body problem can be rephrased as a variational determination of the two- 
body reduced density matrix, subject to a set of A^-representability constraints. The mathematical 
problem has the form of a semidefinite program. We adapt a standard primal-dual interior point 
algorithm in order to exploit the specific structure of the physical problem. In particular the 
matrix-vector product can be calculated very efficiently. We have applied the proposed algorithm 
to a pairing-type Hamiltonian and studied the computational aspects of the method. The standard 
A^-representability conditions perform very well for this problem. 
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I. INTRODUCTION 



It was realized in the 1950's [H |2] that the energy of a quantum many-body system can 
be expressed in terms of the two-body reduced density matrix (2DM), when only one- and 
two-body interactions are present. This insight led to the idea of variationally determining 
the 2DM by minimizing the energy, henceforth referred to as the v2DM method. Once 
the 2DM is known, all other physical properties that can be expressed as one- or two-body 
operators can be extracted. In this way the 2DM effectively replaces the wave function 
and we have "quantum mechanics without wave functions" [3j. Early attempts, however, 
produced unrealistic results [4] and it was soon realized |5j that non-trivial constraints are 
needed to ensure that the 2DM is derivable from a physical wave function. These constraints 
were called iV-representability conditions by Coleman [6] , and Garrod and Percus |7] derived 
two such conditions, the so-called Q and G conditions, which can be expressed as matrix- 
positivity constraints. With these constraints there were some attempts, some of which quite 
successful, to solve this problem numerically in the 1970s [HHTT]. However the method was 
soon abandoned because of the computational cost. Interest in the subject was renewed at 
the beginning of this century, when first Nakata [12] and then Mazziotti [13] realized that the 
v2DM problem can be formulated as a semidefinite program (SDP) for which general-purpose 
primal-dual SDP solvers can be used [H], and they calculated the ground-state properties 
of small atoms and molecules. Primal-dual interior point methods are the "Rolls Royce" of 
SDP algorithms, having several appealing features, but they require a lot of storage and are 
computationally expensive. These early calculations were therefore limited to small systems 
(minimal basis set). Mazziotti |15] then developed an algorithm that transforms the SDP 
into a non-linear optimization program solved by a gradient-only method. This reduced the 
cost of the storage and the basic floating point operations, but at the cost of these nice 
convergence properties of the interior point methods. In this paper we adapt a standard 
primal-dual interior point algorithm [TB] to the specific case of v2DM, in an attempt to 
retain the nice convergence properties, while reducing the storage and computational cost. 
In Sec. [TT] we present an introduction to the theory of A^-representability, v2DM and some 



mathematical properties of the constraints. In Sec. |III| we discuss the representation of the 
problem as a primal-dual semidefinite program, and introduce the method we use to solve 
it. Then we apply the algorithm to a BCS (Bardeen-Cooper-Shrieffer) [T^ or pairing-type 



Hamiltonian in Sec. IV and present the physical results and computational aspects. A 
summary is provided in Sec. |V 



II. VARIATIONAL DENSITY MATRIX DETERMINATION 

When only two-body interactions are present, the Hamiltonian of a physical system can 
be written as: 

H = ^ ta-yO^ay + - ^ Vai3;-ysa'a(i^i3asa^ , (1) 

using second quantized notation where a]^ (oq) creates (annihilates) a fermion in a single- 
particle (sp) state a [IB]. The expectation value of the energy in an arbitrary A^-particle 
state 1^^) can be expressed in terms of the 2DM only, 

E{T) = TtTH^'^= J2 ^^P-nsH^'L. (2) 

a</3;7<5 

with the 2DM defined as: 

r„ft75 = (^^lat^^a^a^l^^) , (3) 
and the reduced two-particle Hamiltonian, 

The idea of v2DM is to determine the ground-state energy and other two- or one-body 
properties by minimizing the energy ^ using the 2DM as a variable. The 2DM is a much 
more compact object than the wave function because one keeps the dimension of two-particle 
(tp) space, no matter how many particles are involved. The problem is that there is no 
straightforward way to know whether an arbitrary matrix in tp-space F is derivable from a 
physical wave function as in Eq. Actually, it is sufficient that T is derivable from an 
ensemble of A^-particle wave functions, and this is called the A^-representability problem [6] . 
Some obvious necessary A^-representability constraints are apparent from the definition (|3]): 

trace condition Tr F 

antisymmetry Ta^.^s 
Hermiticity Fq,^;^^ 

but it turns out that there are many non-trivial constraints needed to ensure that a 2DM is 
physical. 
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A. A^-representability 



The necessary and sufficient conditions for A^-representability are formally known [T^. A 
tp-matrix is A^-representable if and only if, for every two-body Hamiltonian H^, the following 
inequality is satisfied: 

Tr i/(2)r > E^^H,) , (8) 

where Eq{H^) is the exact A^-particle ground-state energy corresponding to the Hamiltonian. 
This is hardly a practical approach, as one needs to know the ground-state energy of every 
two-body Hamiltonian. Therefore one resorts to certain classes of Hamiltonians for which 
a lower bound to the ground-state energy is known. A Hamiltonian class that is used as 
necessary constraint is 

{^^\B^B\^^) > , (9) 

which leads to positivity conditions of linear matrix maps of the 2DM. If we want ^ to 
be restricted to tp-space there are three possible forms of the operator B"^ , leading to three 
conditions on the density matrix: 

a. B"^ = J2af3Po'i3(^a(^^i3 leads to the trivial P-condition: 

V{T) = r ^ , (10) 

which imposes positive semidefiniteness on the 2DM. 

b. B^ = Z ai3 Qai^aaap leads to the Q-condition: 

Q(r) y , (11) 

where the linear matrix map Q is defined as 

f 



N{N-1) 



with 

1 



^ ra/3;7/9 ) (13) 



/3 

the one-body reduced density matrix (IDM), and with 

^ = J2^aP;ap, (14) 

a/3 

the unrestricted trace of the 2DM. 
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c. = J2ai3 QaiiCi^a'^p which leads to the ^-condition: 

g{v) y , (15) 

with the hnear matrix map Q defined as 

= SpsPa-y — '^a5;jf3 ■ (16) 

Another Hamiltonian class for which a lower bound to the ground-state energy is known 
gives rise to the so-called three-index conditions: 

>0 . (17) 



In this article we will use two conditions that come from Eq. (17). 
d. B"^ = J2ai3-f tais-yO-aO-^pO^l leads to the 7i-condition: 

ri(r) y , (18) 

with the linear matrix map 71 defined as 



r 



N{N -I) 

+ (S^C^iss — Si3,^5ys) Pae — {^-yC^aS — ^aC^-fS) Pe/3 + {^K^ai — SaC^lSs) P"/t 
— {SpS^-re — Spt^^/S) PaC + {^-yt^aS — ^ae^-yi) PK " {^Pe^aS — ^at^lSs) P7C 

07; 

+6.ysTai3-e(; — Sj3sTa-y-e(; + SaS^I3r,eC ■ (19) 

e. B'^ = J^ap-y tai3-yCi^a'^\cb-f Isads to the 72-condition 

T2iT) h , (20) 
with the linear matrix map T2 defined as 

T2{T)ai3r,Se(: = (^^|ala^a^4a,a5 + aj^a^ao^aja^l^^) (21) 

= {^aS^Pe — ^ae^ps) P-yC + SyC^al3;Se — SaS^yt;Cl3 + Sl35^yt;(a + Sae^yS;CP ~ ^Pt^yS;Ca ■ 



The optimization problem that we have to solve can be summarized as: 

mm Tr TH^^^ , (22) 

under the condition that 

Tr r = ^ , (23) 

£(r)^0 V£g {P, 2,6^,7^, Ts} . (24) 

B. Hermitian adjoint maps 

For the following it is useful to introduce the Hermitian adjoints of matrix maps intro- 
duced in the previous section. The Hermitian adjoint maps are defined through: 

Tr A(r)A = Tr C\{A)V , (25) 

in which A is a matrix of the same dimension as the image of the map in question i^e.g. 
a three-particle matrix for a T\ map, eic), and the traces sum over the appropriate indices. 
The V and Q, maps are Hermitian, so they are identical to their Hermitian adjoints. For 



the other maps however this is not the case. Using Eq. (25) the Hermitian adjoint of the Q 
map can be shown to have the form: 

(^)a/3;7<5 = - \ \^fi^-^'^1 ~ ^aS^p.^ - dfi^AaS + 5^7^/35] (26) 

in which a particle-hole matrix A is mapped onto tp-matrix space and 

A-a-f = ^ ^ ^Q!A;7A • (27) 
A 

The 7i-operator maps a tp-matrix onto a three-particle matrix, so its Hermitian adjoint has 



to map a three-particle matrix A onto tp-space. Solving Eq. (25) with £ = 7i one finds 
that: 

2 

T'l (^)a/3i75 = {^a^ihs " ^cS^p^) Tr A + A^p,^s (28) 



2(iV- 1) 
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with 



In the same way one can derive ior C — % that 



2{N-1) 

Asa;l3'y ~ Agp-Q.^ — A^a;/3S + A^p-^S 



(29) 
(30) 

(31) 



with this time A a matrix on two-particle-one-hole space and 



Anft-^ti = ^ ^ A 



Q:/3;7d — ^afiX\^5X 
X 



(32) 
(33) 
(34) 



III. PRIMAL-DUAL SEMIDEFINITE PROGRAM 

The variational method described in the previous section can be formulated as a primal- 
dual semidefinite program. A general 2DM, describing an A'"-particle system can be expanded 
in an arbitrary orthogonal basis {/'} of traceless matrix space as 



N{N -I) 



(35) 



M(M- !)■ 

with M the dimension of single-particle (sp) space, and the unit matrix on tp space defined 



as 



The energy of the system can be written as a function of the 7's as 



M M - 1 ^ J 



(36) 



(37) 



M(M-l) 

Because the necessary A'"-representability conditions can be written as linear homogeneous 
matrix maps of F, we can also write them as a function of the 7's: 

N{N - 1) 



^(r) = ]g^£(i.p) + E^.^(/') to- 



ps) 



If we now consider the direct sum of the hnear spaces associated with the maps and define 
the block matrices: 

N{N-1) 



(Itp) and = if) , (39) 



M(M - ^, 

^ ' k k 



then we can formulate v2DM as a standard dual semidefinite program 

min 7^/i on condition that Z = + 7jM* ^ , (40) 



in which = Tr H^"^^ f\ The primal problem corresponding to (40) optimizes the matrix- 
variable X, the problem being defined as: 

max (-Tr Xu°) on condition that Tr Xu' = h' and X ^ . (41) 

X will be a block matrix because the w-matrices are block matrices. The primal-dual gap 
rj is defined as the difference between the primal and the dual cost function for a certain 
primal-dual point (X, Z): 

= + Tr u^X = ^ 7iTr Xu' + Tr Xm° = Tr > , (42) 

i 

as X and Z are positive semidefinite matrices. We can see that the smallest value of rj 
will be reached when both the primal and the dual problem are optimal. It can be proven 
that if the primal and the dual problem are both strictly feasible, then the primal-dual 
gap vanishes at their solution pH]. This means that the primal-dual gap can be used as a 
convergence criterion for the algorithm. Even better, at any point during the optimization, 
the error on the current value is limited from above by the primal-dual gap. Note that in 
our previous implementation f^U], a dual-only algorithm was used. The properties of the 
present primal-dual method can lead to a serious reduction in computation time since we 
can stop the algorithm at a prescribed error estimate. 

A. Equations of motion 

There are several known methods to solve a semidefinite program. In this paper a path- 
following interior point method is used. The central path is defined as the set of primal-dual 
points for which 

XZ = ^Is^p , (43) 
n 
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with n the total dimension of the X and Z matrices and Isup the direct sum of the unity 
matrices on the different constraint spaces: 

Isup = Ifc . (44) 

k 

In the path- following algorithm [TB] we try to follow the central path, reducing the primal- 
dual gap along the way. Consider a primal-dual point (X, Z) on the central path with 
primal-dual gap rj. We want to know what is the primal-dual point on the central path with 
primal-dual gap scaled down with a factor p. Rephrasing, we are looking for the (Ax, A^) 
that solve: 

(X + Ax)(^ + Az) = ^l,,p. (45) 

There are several ways to symmetrize these equations. Using the method proposed by [T6] . 
two equivalent equations (called the dual and the primal) are obtained, i.e. one has to solve 
the equations 

(dual) : Ax + D-^AzD-^ = ^Z'^ - X , (46) 

n 

(primal) : Az + D Ax D = ^X'^ - Z , (47) 
and under the condition that: 

Tr Axu' = and A^ = ^ (^7)^^' , (48) 

i 

and with 

D{X,Z) = X-^ (^X-^ZX'^y X--2 . (49) 
1. Solution to the dual equation 

In order to obtain the primal-dual direction (Ax, A^) , the dual equation (46) is first 



projected onto the space spanned by the non-orthogonal basis (which we will call U- 



space). With B denoting the right-hand side of (46) and making use of Eq. (48) we obtain 



V (Tr D-\W-\') A7j = Tr Bu' , (50) 



which can be seen to be a symmetrical, positive-definite linear system and as such can be 
solved iteratively using the linear conjugate gradient method. This can be done without 



explicit construction of the dual Hessian matrix T-i^ or any reference to the non-orthogonal 
basis set {u^}- This is because can be seen as a map from traceless tp-matrix space onto 
itself, by using the Hermitian adjoints of the linear maps C Consider an arbitrary traceless 
tp-matrix: 

6 = ^6,/^-. (51) 



Using (25) and the fact that the £'s are linear and homogeneous we obtain that the image 



of e under the dual Hessian map can be written as: 



Tr 



(52) 



in which the Dk are the blocks of the D matrix corresponding to the different constraints 
Ck, and Ptt stands for the projection operator onto traceless tp-matrix space: 

2Tr A 



M(M- 1) 



1 



tp 



(53) 



2. Solution to the primal equation 



The solution of the primal equation (47) is obtained in the same manner, by projecting 



this equation onto C-space, the orthogonal complement of W-space. With B denoting the 



right-hand side of the equation (47) and making use of Eq. (48) one gets 



V (Tr D D d) 5xj = Tr Be' 



(54) 



«5 



where we have used 



A 



X 



(55) 



This is again a symmetrical positive-definite system of linear equations that can be solved 
iteratively using the linear conjugate gradient method. As with the dual equation it can 
be solved without explicit construction of the Hessian matrix ^ or any reference to the 
basisset {c*}, because 'K^ can be seen as a map from C-space onto itself. For an arbitrary 
matrix in C-space: 

e = 5^e,c\ (56) 
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the image of e under the primal Hessian map is 



n^e = Pc [DeD] 



(57) 



in which Pc is the projection onto C-space. This projection can be executed quickly by 
using the inverse of the overlap matrix of the W-space basis vectors. Suppose we have an 
arbitrary block matrix A of the same dimension as X and Z. First we project it onto the 
space spanned by the basis {u^,u^} = The projected matrix A' reads as: 



(58) 



where the overlap matrix S appears because of the non-orthogonality of the basis. Due to 
the special properties of the linear matrix maps C that determine the basis matrices u", 
the inverse overlap matrix can also be considered as a map from tp space onto itself (see 
Appendix |A] and [B] for the actual analytic expression of this map). The projected matrix A' 
can now be written in block-matrix form as: 



(59) 



To project A onto W-space we still have to remove the component along the u°-matrix: 

'Tr u^A' 



PuA = A' 



Tr v?u^ 



u 



(60) 



Since C-space is the orthogonal complement of the W-space, the desired projection of A onto 
the C-space is simply given by 

PcA = A- PuA . (61) 



B. Outline of the algorithm 

In this section a short outline of the algorithm will be presented. The first step is to 
initialize the primal-dual variables, after which they are directed towards the central path. 
Then the actual minimization of the primal-dual gap takes place, which is done in a predictor- 
corrector loop. 
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1. Initialization 



We need a feasible primal-dual starting point. An initial feasible dual point Z^'^\ i.e. a 



matrix that satisfies the inequality (40), is easily found by setting 



(62) 



which corresponds to setting al the 7j's equal to zero. A feasible primal starting point will 



have to satisfy Eq. (41). To construct such a point we take a completely random matrix X 



and project it onto a matrix X' for which 



Tr X'u' = h' 



(63) 



This is again achieved using the inverse overlap matrix of the {u""} basis, 



X' = X - J] (Tr Xu" - /i") S-^u^ 



The last term on the right-hand side can be computed as: 



(64) 



(65) 



At this point, X' satifies the equality constraint (63), and one just has to add vP , with a 



positive scaling factor that is large enough to ensure positive semidefiniteness: 



= X' -t- ^ . 



(66) 



2. Centering run 



Before the actual program can be started, a couple of centering steps have to be taken. 



which is done by solving the equations (46) and (47) with u = 1. The purpose is to go 



sufficiently near the central path, without bothering about the primal-dual gap. In a first 



step, Eq. (50) which has the smallest dimension, is solved using the conjugate gradient 



method, and the dual solution Az is obtained. The primal solution Ax then follows from 



the dual equation (46) by substitution. For these initial centering steps, both linear systems 
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are so well conditioned that hardly any iterations are needed for convergence. As a measure 
for the distance from the center we use the potential [r4j : 



= -IndetX-lndetZ , (67) 
which is minimal (for points with the same primal-dual gap rj = Tr XZ) on the central path 



for which Eq. (43) is satisfied: 



$(X^Z") = -nln^ . (68) 



When the potential difference (which is always positive) 



*(X, Z) = ^{X, Z) - <I>(X^ Z^) 

= n In Tr XZ — n In n — In det X — In det Z , (69) 



is sufficiently small, the centering run is stopped. 



3. Predictor- corrector run 



In this part of the program the primal-dual gap is minimized by alternating predictor and 
corrector steps. A predictor step tries to reduce the primal-dual gap by solving the equations 



(46) and (47) with z/ = 0. This is done in exactly the same way as for the centering run, 



by first solving (50) for Az, then substituting into (46) to obtain an approximate primal 
step Ax. The final primal step is obtained by solving (54) using the conjugate gradient 
method with the approximate Ax as a starting point. Note that when the primal-dual gap 
decreases, the condition number of the primal and dual Hessian matrices increases and more 
iterations are needed before convergence is reached. One can adjust the convergence criteria 
of the primal and dual conjugate gradient loops, in order to minimize the combined number 
of iterations. 

At this point we have a predictor direction (Ax, A^). The logarithmic potential 0(a) = 



\1/(X -|- aAx, Z + aAz) in the predictor direction (see Eq. (69)) can be simply evaluated for 
any value of a by precomputing the eigenvalues Xf of X~5AxX~^ and Af of Z^^ AzZ~K 
One then has 



[a 



^(X, Z) + In [1 + a{cx + cz)] - ^ ln(l + aXf) - ^ ln(l + aXf) , (70) 
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FIG. 1. The ground-state energy as calculated by v2DM(PQG) and the Richardson-Gaudin 
equations (RG), together with the pair occupation in the groundstate by v2DM(PQG), as a function 
of the pairing interaction strength g. 

where 

cz = -Tr XAz and cx = -Tr ZAx , (71) 

T] T] 

With a standard bisection method one can now compute the stepsize a corresponding to 
the maximal deviation from the central path we want to allow. 

After the predictor step, a corrector step is taken, which is equivalent to the centering 



step described previously (see Sec. IIIB2). The alternation of predictor and corrector steps 



continues until the primal-dual gap is smaller then the desired value. 
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FIG. 2. The difference between the ground-state energy calculated by v2DM with various con- 
straints, and the exact solution, as a function of pairing strength g. 

IV. APPLICATION TO THE BCS HAMILTONIAN 
A. The BCS Hamiltonian 



The algorithm introduced in Sec. Ill is applied to the BCS Hamiltionian [T^. The BCS 
Hamiltonian is an interesting system that models the competition between a single-particle 
operator and a schematic pairing interaction: 

H = Y^ eia^tti^ -gY^ al^al^ttj^aj^ . (72) 

icr ij 

Here the single-particle levels are denoted with an index i = 1, . . . , M, and the up (down) 
spin as a =t (4-)- When the pairing strength g is small compared to the single-particle level 
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spacing, the energy is minimized by filling up the single-particle orbitals up to the fermi level. 
With increasing g however, it becomes advantageous to form pairs, i.e. it is energetically 
favorable to maximize the ground-state occupation of the fermion pair state ^2 - al^al^- This 
problem is hard to solve using standard perturbative methods as these tend to break down 
when pairs are formed. An exact solution based on the Bethe-ansatz exists for this problem, 
however, and involves solving a system of non-linear equations pT]. These equations are 
notoriously difficult to solve because, for certain critical values of g, the equations become 
singular. Several approaches have been suggested for solving these equations [221 123] • In 
this paper we follow the approach recently proposed by De Baerdemacker [21]. The exact 
ground-state energies as a function of g are compared to the v2DM results calculated within 
the present formalism. 



B. Results 



We have studied the Hamiltonian Eq. (72) with M = 12 doubly degenerate equidistant 



single-particle levels and = 12 fermions, and g ranging from to 5 in steps of 0.01. v2DM 
calculations were performed with respectively PQG, PQGTi and PQGT1T2 constraints. The 
resulting ground-state energy is compared to the exact solution in Fig. [1} For all values of g 
the agreement is already remarkably good at the PQG level. To appreciate how the result 
improves when constraints are added the difference between the various v2DM results and 
the exact solution is plotted in Fig.|2} Note that the difference is always negative, since v2DM 
provides a variational lower bound. As one observes, all approximations describe exactly 
the non- interacting small- (7 limit. When g becomes larger, there is competition between 
different types of ground states and the performance of PQG gets worse up to (7 ~ 1.4. For 
larger g the PQG result becomes better again. In fact, we checked (by omitting the single- 
particle piece) that also the g ^ 00 limit becomes exact for PQG, which is a peculiarity of 
the schematic pairing force. The PQGTi results show that the Ti constraint only becomes 
active around g = 2.5, and ensures faster convergence to the exact g ^ 00 limit. Somewhat 
surprisingly, adding the T2 condition is sufficient for obtaining the exact solution at all values 
oi g. 
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C. Computational Performance 

Some of the computational aspects of the algorithm are worth pointing out. It is inter- 
esting to see e.g. that depending on the pairing interaction parameter g, the convergence 
properties of the algorithm change. In Fig. |3] the joint number of predictor and corrector 
steps needed for convergence, is plotted as a function of g. One observes a sharp peak at 
fairly small g, just when the perturbative regime is left and the structure of the ground state 
changes. For g = 0.25, which is at the position of the peak in Fig. [3| we have plotted in Fig. 
|4]the number of conjugate gradient iterations needed for convergence, of both the dual and 
the primal linear system, as a function of the primal-dual gap 77. As expected, the number 
of iterations for the dual problem increases with decreasing primal-dual gap, as the linear 
system grows ill-conditioned. The primal conjugate gradient loop only becomes active for 
small values of 77. This signals that the numerical stability becomes too small to generate 
a high quality approximation for using the Az obtained in the dual conjugate gradient 
loop. Anyway, the needed number of primal iteration remains insignificant compared to the 
dual ones, for all values of rj. The situation a.t g = 0.25 is the worst case. For larger values of 
g, where the number of predictor-corrector steps is smaller and approximately constant (see 
Fig. [3]), the number of conjugate gradient iterations is also drastically reduced. A typical 
behaviour is plotted in Fig. [5] for (7 = 4. 



V. SUMMARY AND DISCUSSION 

Interacting quantum many-particle systems lie at the heart of most issues in condensed 
matter, molecular/atomic and nuclear physics. Their analysis may be rephrased as the 
problem of minimizing the energy, expressed as a linear function of a two-body density 
matrix, subject to the A^-representability constraint that the 2DM can be derived from a 
physical A^-particle system. By working solely with the 2DM, rather than with the A^- 
particle wave function itself, the problem of the exponentially exploding dimension of A^- 
particle Hilbert space with increasing A^ is circumvented. The complexity of the problem is 
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FIG. 3. Number of predictor and corrector steps needed for convergence, as a function of the 
pairing strength g in v2DM(PQG). 

shifted, however, to the characterization of the iV-particle representable 2DM's. In practice, 
a limited set of necessary but not sufficient conditions for A^-representability are imposed 
during the minimization, resulting in a strict lower bound to the energy, which converges to 
the exact energy when more and more A^-representability conditions are imposed. 

Commonly used A^-representability conditions impose the positive semidefiniteness of a 
set of linear matrix functionals of the 2DM. In this way the quantum many-body problem 
is converted into a well established field of optimization techniques called semidefinite pro- 
gramming. Standard packages for SDP, however, fail to take into account properties of the 
physical problem that can be exploited. 

Using specific mathematical properties of the constraints for the v2DM problem, we 
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FIG. 4. Number of primal and dual conjugate gradient iterations needed for convergence, as a 
function of the primal-dual gap r] for g = 0.25 in v2DM(PQG). 



have adapted a standard primal-dual interior point method to be computationally cheaper, 
both in storage as in floating point operations. We make extensize use of the algebra 
of linear matrix maps to calculate efficiently some intermediate quantities. During the 
Newton minimization procedure, a new direction in 2DM space is found iteratively using the 
conjugate gradient algorithm, thereby exploiting the fact that the product of the Hessian 
with a 2DM is considerably cheaper for the physical problem at hand than in a general 
situation. 

As an example we have applied the algorithm to a BCS-type Hamiltonian. We found 
that the standard constraints work very well for this kind of problem. The computational 
performance of the method was analyzed, and it was shown that the convergence behaviour 
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FIG. 5. Number of primal and dual conjugate gradient iterations needed for convergence, as a 
function of the primal-dual gap r] for g = 4 in v2DM(PQG). 

is dependent on the value of the pairing strength parameter. As in our previous algorithm 
[20] the method slows down near the solution, because the matrices involved become ill 
conditioned. The present primal-dual algorithm allows to control this since the primal-dual 
gap provides an upper bound to the remaining error. Therefore the algorithm can be stopped 
when the required accuracy is reached, saving many unnecessary iterations. 
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Appendix A: Calculation of the overlap-matrix map 



The overlap matrix of the non-orthogonal basisset is defined as: 

Sai3 = Tr . 

Using the Hermitian adjoints of the linear maps C we can rewrite this as: 

5«;, = ^Tr \cl{CUn)f 



(Al) 



(A2) 



in which {f^} is an orthogonal basis of tp-matrix space. This means that the overlap matrix 
can be seen as a linear map from tp-space onto itself, whose action onto a tp-matrix F is: 



'5(r) = $^4(A(r)) . 



(A3) 



It turns out that this map can be written as a generalized Q map, which is defined as: 

Q{a, b, c) (X)ai3-^S ^ (^^al3;-f5 + b {Sa-ySps — ^aS^jB-i) T — C (^a^r^^ — 5 p^V — ^aS^ Py + ^pS^ay) ■ 

(A4) 



This is like a Q-map (12) but with general coefficients (a, 6, c). The proof is somewhat 
tedious and proceeds by considering every Ck separately. 



1. 



It is trivial to see that p2(r) = r and that this is a generalized Q map with coefficients 



a=l 6 = c = 0. 



(A5) 



2. Q2 



To reexpress we first calculate the various pieces. 



M-N-l 



_ N{N - 1) 
(M- A)(M-A^-1) 



N- 1 
f . 



A^(A- 1) 

Substitute into Eq. (12) leads once again to a generalized Q map with coefficients: 

4A2 + 2N- ANM + - M 



a = 1 



N^{N-iy 



2N -M 
(A - 1)2 



(A6) 
(A7) 

(A8) 
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3. g^g 



With the same strategy one finds on the basis of Eq. ( 26 ) and 



a7 } 



that substituting into (26) leads to another generahzed Q map with coefficients: 

2A^ - M - 2 

a = 4 6 = C-- 



(A9) 



(AlO) 



4. VTi 

The needed terms are now: 



M -N-2 



Ti (r)„, 
ti(r) 



A^- 1 

(M- A^-2)(M- A^- 1 



_ A^(A^- 1) _ 

(M-3)(M-2A^) 



A^(A^-l) 
(M - 2)(M(M - 1) - 3A^(M - A^)) 
Ar(A^- 1) 



A^- 1 



r 



(All) 
(A12) 
Fq,^ , (Al3) 
(A14) 



and substitution into Eq. (28) leads to the coefficients 



a = M - 4 , 

M3 - 6M^N - 3M2 + 12MAr2 + 12MAr + 2M - ISA^^ - QN^ 
M2 + 2A^2 _ 4j^^jY _ jvf + 8A^ - 4 



2(A^- 1)2 



5. T^r2 



Finally, needed for the calculation of the last map are: 

f 



7-2 (F), 



M-N- 

^2 (r)„^.^5 = ^ _ ]_ + ^pS^a-i - (M - 2)F„5.^/3 , 

'M(M - A^) - (A^ - 1)(M - 2)' 



r2(F) 



ay 



N -I 



(A15) 
(A16) 
(A17) 

(A18) 
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which, when substituted into Eq. (31) gives the following coefficients: 

2 2A^2^ (M-2)(4A^-3) 



5M-8 



M2 



(A19) 



2(iV-l)2 

The overlap-matrix map is just the sum of the various terms obtained, and hence also a 
generalized Q map with rather complex coefficients. 



Appendix B: Inverse of generalized Q map 



The inverse of a generalized Q map can be shown to be another generalized Q map. 
Consider for brevity the notation: 



Q{a,b,c)iT)=Q 



then applying partial trace operations on Eq. (A4) leads to: 



a + M{M - l)b - 2{M - l)c 



07 



a - c{M - 2) 



6(M-1) -c 



a + M{M - l)b - 2{M - l)c 



Upon substitution into Eq. (A4) and solving for T one obtains. 



T = Q-\a,b,c)iQ) = Qia',b',c')iQ) 



where 



b' 



ba + bcM - 2c2 



a [c(M - 2) - a] [a + bM{M - 1) - 2c(M - 1)] 
c 



(Bl) 



(B2) 
(B3) 

(B4) 

(B5) 
(B6) 
(B7) 



a [c(M - 2) - a] ' 

These are important relations since they allow to evaluate the action of the inverse overlap 
matrix on a tp matrix as fast as a Q map. i.e. at a computational cost which is negligible 
compared to the other matrix manipulations. 
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